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Abstract — An  overview  of  the  multipath  expansion  method  of 
solving  the  Helmholtz  wave  equation  to  describe  the  underwater 
sound  field  for  a  fixed  point  source  in  a  plane  multilayered  medium 
is  presented.  The  approach  is  then  extended  to  account  for  hori¬ 
zontal  variations  in  bottom  depth,  bottom  type,  and  sound  speed  in 
the  stationary  phase  approximation.  Comparisons  of  model  results 
to  a  limited  number  of  measured  data  sets  and  standard  propaga¬ 
tion  codes  are  presented. 

Index  Terms — Acoustics,  multipath  expansion,  stationary  phase, 
underwater  sound  propagation,  wave  equation. 

I.  Introduction 

HIS  paper  presents  an  overview  of  the  derivation  of  the 
propagation  loss  model  of  Leibiger  [1],  [2],  and  then  ex¬ 
tends  this  result  to  shallow  water  environments  by  accounting 
for  changes  in  bottom  depth,  bottom  type,  and  sound  speed  with 
range.  Comparisons  between  calculated  and  measured  results 
are  presented. 

A  derivation  based  in  part  on  the  work  of  Leibiger  was 
published  by  Weinberg  in  1975  [3],  [4],  Recent  developments 
by  Weinberg  extend  the  use  of  ray  methods  through  the  use 
of  Gaussian  ray  bundles  [5]  but  do  not  utilize  the  integration 
method  discussed  in  this  paper.  A  comparison  of  results  for 
a  transmission  loss  (TL)  calculation  at  3500  Hz  using  the 
Integrated  Mode  method  described  herein  and  Weinberg’s 
Navy  standard  Gaussian  RAy  Bundling  (GRAB)  model  is  also 
presented. 

The  calculation  of  sound  propagation  in  the  low  frequency 
portion  of  the  sonic  band  (20  to  200  Hz)  has  been  addressed 
using  ray  theory,  parabolic  equation  (PE)  methods,  and  normal 
mode  codes,  with  which  this  approach  agrees.  However,  the  cal¬ 
culation  method  described  herein  is  also  viable  at  higher  fre¬ 
quencies,  from  200  Hz  to  above  50  kHz.  The  need  for  calcula¬ 
tions  which  are  continuous  across  the  operational  band  of  Navy 
platforms  and  systems  in  shallow  water  provides  motivation  for 
the  current  work. 

An  overview  of  the  Leibiger  derivation  is  presented  in 
Section  III,  followed  by  an  extension  to  range-varying  media 
in  Section  IV.  Section  V  provides  a  comparison  of  results  to  a 
limited  number  of  measured  data  sets  and  other  propagation 
codes. 


Manuscript  received  April  15.  2003;  revised  June  15.  2004;  accepted  July  15, 
2004.  Associate  Editor:  W.  Carey. 

The  author  is  with  the  Naval  Undersea  Warfare  Center,  Newport,  RI  02841 
USA. 

Digital  Object  Identifier  10.1 109/JOE.2004.834167 


II.  Preliminaries 

We  define  sound  speed  c  to  take  values  within  the  vertical 
plane  of  the  ocean  defined  by  the  positions  of  the  source 
and  receiver  of  the  acoustic  energy.  We  then  define  the  total 
wavenumber  k  to  be  equal  to  oj/c,  where  u  is  the  radian 
frequency  of  the  source  and  denote  the  horizontal  and  vertical 
vector  components  of  the  wavenumber  k  by  kr  and  kz,  respec¬ 
tively,  where  r  denotes  range  and  z  denotes  depth.  We  define 
6  to  be  the  angle  of  the  total  wavenumber  with  the  horizontal, 
and  denote  the  vertical  phase  function  £  as 


Assuming  a  region  which  is  homogeneous  in  range,  we  evaluate 
£  and  its  derivatives  by  partitioning  the  local  sound  speed  pro¬ 
file  (SSP)  into  linear  segments  in  the  usual  way,  thus  defining  a 
horizontally  stratified  ocean.  In  horizontally  variable  media,  we 
likewise  partition  the  range  of  interest  into  segments  for  which 
1)  sound  speed  variation  with  depth;  2)  bottom  slope;  and  3) 
bottom  type  are  considered  invariant.  Thus,  calculations  are  al¬ 
ways  performed  within  a  rectangle  in  which  the  sound  speed 
variation  is  the  same  linear  function  of  depth  throughout  the 
region. 

The  horizontal  range  traversed  in  propagating  through  a  layer 
from  z\  to  Z‘2  with  constant  sound  speed  gradient  g  is  given  by 

Ar  =  — — Asin0  (2) 

9 

where  a„  denotes  the  vertex  sound  speed  at  which  all  propaga¬ 
tion  is  momentarily  horizontal  (6(z)  =  0). 

Substituting  kz  =  \J (ut/c)2  —  k%  into  (1),  using  the  fact  that 
kr  =  (u/cv),  performing  a  change  of  variable  dz  =  (1  / g)dc, 
and  integrating  yields 


Thus  the  change  in  vertical  phase  in  traversing  a  horizontal  layer 
from  Z]  to  Z2  is  computed  as 

A£  =  £(c(z2)  -  £(c(zi)).  (4) 
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N  ote  that 


_d_ 

dkr 


kz(z,kr ) 


and  so 


differentiating  (3)  with  respect  to  kr  yields 


For  k2z  <  0  we  maintain  positive  values  for  phase  by  adopting 
the  convention  kz  =  — t|Az|,£  =  i|£|.  Substitution  of  the 
leading  term  of  £  from  (8)  into  the  homogeneous  form  of  the 
^-dependent  wave  (7),  letting  Z  =  (z  -  a)1/2y,  and  changing 
the  independent  variable  to  £  yields  Bessel’s  equation  of  order 
1/3  [10] 


<7£  1/ /  1  \  dkz  u j\  1  ckz 

dkr  9  \  \  u>/c+  kz )  dkr  kr  J  g  kr 


dry  1  dy 
d?+~^  + 


(9) 


M  ulti plying  and  dividing  by  e,,  and  utilizing  Snell’s  Law  yields 

d£  _  _  _^sin0 

dkr  g  c„  kr  g 

For  use  in  the  sequel,  we  note  that 


_d_ 

dkr 


k^dz  = 


dk , 
<7A, 


dz  = 


'21 


22  k 

—dz. 
kz 


(5) 


III.  Overview  of  the  Leibiger  Derivation 

Wenow  give  an  overview  of  the  Leibiger  derivation  in  prepa¬ 
ration  for  extending  to  horizontally  variable  media;  a  complete 
discussion  of  the  original  work  may  befound  in  [1],  Derivations 
of  the  acoustic  wave  equation  and  further  discussion  of  various 
solutions  can  also  be  found  in  multiple  references,  see,  for  ex¬ 
ample  [6]— [9]. 

Assuming  sound  pressure  ip  has  no  azimuthal  dependence, 
the  wave  equation  can  be  expressed  in  cylindrical  coordinates 
as 


d2ip  1  dip  d2ip  _  1  d2ip  .  . 

dr 2  r  dr  dz2  c2  <9i2 

where  r  represents  horizontal  range,  and  depth  z  is  increasing 
down  from  the  ocean  surface.  After  a  standard  separation  of 
variables,  the  depth  dependent  function  Z(z)  is  given  by 

Z"(z)  +  k2z(z)Z(z)  =  S(z  -  zq)  (7) 

where  kz(z)  is  the  vertical  component  of  the  local  wavenumber, 
and  weassumea  harmonic  point  sourceof  unit  i  ntensity  at  depth 
zq.  Weconstruct  solutions  of  the  depth  dependent  wave  equation 
from  the  Green’s  function  for  the  associated  boundary  value 
problem  (see,  for  example,  [11],  [12]).  At  this  point,  we  drop 
the  subscript?’  from  the  horizontal  wavenumber  to  simplify  the 
notation,  i.e.,  in  the  sequel,  the  unsubscripted  variable  A  implies 
kr.  Let  z  =  a  be  a  turning  point  for  the  wavenumber  k,  thus 
k2(a,  k )  =  0.  We  assume  that  k2z  can  be  expressed  as  a  linear 
function  of  z  in  a  neighborhood  of  a,  and  let 

kz(z,  k )  =  M2(z  —  a)  +  o(z  —  a ) 


Thus  cylinder  functions  of  order  1/3  provide  solutions  for  Z, 
for  example, 


Z  =  (z-a)1/2J1/3(0.  (10) 

Noting  that  z  -  a  is  proportional  to  /At1,  a  set  of  allowable 
solutions  for  propagation  toward  and  away  from  a  when  com¬ 
pared  with  time  dependence  eiujt  are  given  by 

Z1  =  K1]JJ-H[1/3(0  and 

z2  =  k2]IJ-h[%(0  (11) 

where  k2  >  0  is  assumed  and  7^  and  I<2  are  suitable  constants. 

The  leading  terms  of  the  asymptotic  expansions  of  7T1/3(£) 
and  7T^(£)  when  arg(£)  >  -7 r  are  given  by  [10] 

Hi%(0  ~  5w/12)_  (12) 

Substitution  of  (12)  into  (11)  yields,  with  the  appropriate  choice 
of  constants 


Zx  ~  \[^e-i5*/12^=e*  and 

V  7T  vA, 

Z2  ~  S^e-™l12^=e-*  (13) 

which  are  valid  away  from  the  turning  point  z  =  a  on  the  side 
where  Icos^l  <  1  (they  become  infinite  when  A2  =  0).  The 
Bessel  function  expressions  (11)  allow  the  solution  of  the  dif¬ 
ferential  equation  to  be  approximated  over  a  specific  ^-interval, 
both  near  and  away  from  a,  and  they  agree  with  the  exponential 
solutions  (12)  in  regions  where  they  are  oscillatory. 

To  simplify  the  notation,  we  choose  the  constants  A"i  and  I<2 
so  that  the  depth  dependent  function  is  simply 

Zi  =  ve* 

Z2  =  ve~ii  (14) 


where  M  is  a  positive  constant  and  both  A  and  a  are  assumed 
real.  Integrating  from  a  to  an  arbitrary  depth  z  >  a  to  obtain 
the  elapsed  vertical  phase  £(z,  A)  yields 

£(2,  A)  =  f  kz(z ,  k)dz  =  \m(z  -  cl)3/2  + -  (8) 

J  rv  ^ 


where  v(z, k)  ~  (1  l\fkz)<  and  determine  the  linear  combina¬ 
tion  of  solutions  Z\  and  Z2  which  solves  the  boundary  value 
problem  in  the  ^-domain. 

For  the  purpose  of  this  paper  we  defer  the  discussion  of  rig¬ 
orous  formulations  for  the  boundary  conditions  and  assume  that 
the  reflection  loss  and  phase  change  at  the  upper  )  and  lower 
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(z2)  boundaries  can  be  characterized  as  a  function  of  grazing 
angle  8  by  complex  reflection  coefficients 

R1  =  \R1\ei*1  and 

R2  =  \R2\e^2  (15) 

respectively,  and  thus  that  the  fundamental  wave  functions  ui 
and  u2  satisfy  the  upper  and  lower  boundary  conditions,  respec¬ 
tively,  i.e. 

u\  =  Z\  +  R\Z2 
u2  =  Z2  +  R2Zi. 

Substituting  the  fundamental  solutions  (14),  then,  we  have 

«1  =  v(z,  k)  (V':  r  k>dz  +  Rie+i  r1  k ^  and 

u2  =  v(z,  k )  (e+if~  2  k'-dz  +  R2r;f  k .  (16) 

We  form  the  Green’s  function  solution  from  the  upper  and  lower 
fundamental  solutions 

hJSL ) '  ,  >  °  (17) 

whereto  is  the  depth  of  the  source  and  W(k )  istheWronskian 

W(k)  =  ui(z)u'2(z)  —  u2(z)u[(z).  (18) 

The  Green’s  function  construction  extended  to  the  range  do¬ 
main  yields  the  pressure  <p  as  a  contour  integral  in  the  complex 
k-plane 

<p(r,  z)  =  2  (j>  HQ2\kr)Z(z,  Zo,k)kdk  (19) 

where  the  integral  encompasses  all  of  the  zeros  of  the  Wron- 
skian. 

We  develop  a  practical  expression  for  ip(r,z)  by  formulating 
the  Wronskian  W(k)  from  the  basic  asymptotic  forms  of 
uxiz^k)  and  u2(z,k )  given  in  (16),  which,  after  substituting 
into  (18),  and  performing  some  algebra,  yields 


To  extract  the  roots  of  W0(k),  consider  first  a  propagation 
medium  in  which  boundary  reflections  are  without  loss,  i.e.,  let 
|i?i|  =  \R2\  =  1.  In  this  special  case,  the  zeros  of  W0(k) 
are  given  by  real  values  of  k  for  which  the  total  phase  change 
in  traversing  the  medium  in  depth  from  upper  turning  point  to 
upper  turning  point  is  an  integer  multiple  of  2n,  i.e., 


/  k,dz  —  <f>i  —  <f>2  =  2mr,  n  =  1,2, _  (22) 

The  nonidealized  ocean  propagation  medium  is  lossy,  re¬ 
quiring  |i?ii?2|  <  1-  This  leads  to  nonzero  imaginary  parts  sn 
of  the  eigenvalues  kn,  which  are  found  by  expansion  of  (22) 
about  the  real  part  of  k.  To  expand  kz(z,k2)  around  Re(£r)  to 
the  first  order 


P2  k,dz 

=  -2ieJ^  ~ 


W(k )  =  -2  ie- 


W0(k ) 


2y/(<. u/c)2  -  k2  kz(z,k 2) 


k2(z,k2)  =  kz(z,Re(k2))  -  k  (f^fc2))(fc  -  Re(fc)). 

But  k  -  R e(k)  =  ie,  so  we  have 

1  -  =  0,  or 

l  -  ronSis5)J'  =  o. 

Combining  this  result  with  (5),  the  imaginary  part  of  the  eigen¬ 
value  for  mode  n  is  then  given  by 

_  lnQRiR2|) 

”  Rc(Mkn)) 

as  can  be  verified  by  substitution  of  the  above  into  the  expression 
for  W0(k)  (21). 

A  Ithough  somevariations  of  the  boundary  value  problem  may 
include  propagation  of  an  infinite  number  of  normal  modes  in 
the  ocean  medium,  in  practical  applications  a  possibly  large  but 
finite  number  of  modes  suffices  to  produce  tp(r,  z).  E  valulating 
by  the  calculus  of  residues  yields 


1  ID  |  |  D  I  *($1  +®2-2  fZ2  kztte)  ni. 

Wo(k)  =  1  -  |i?i||i?2|e  J=1  .  (21) 

The  normal  modes,  which  correspond  to  the  zeros  of  W0(k), 
take  sequential  values  A;  =  kn,n  =  1,2, . .  .where?!,  is  the  mode 
number.  In  general,  the  kn  are  complex  with  small  imaginary 
part 

kn  —  Rc( k,,  J  -J-  7cr, ,  En  0,  |c„  |  |Re(/u’„)|. 


z)  =  —in  ^  H^\kr)G(z,  k)k  e 


I  k=He(kn) 

(24) 


where  G(z, k)  =  ( u1(z,k)u2(z,k)/Rc(k )).  Note  that  all  ex¬ 
pressions  in  the  above  sum  may  be  evaluated  using  the  real  part 
of  the  eigenvalue  with  the  exception  of  the  range  dependent 
Hankel  function  in  which  the  imaginary  en  provides  an  expo¬ 
nential  attenuation  factor  e£nr  which  accounts  for  the  boundary 
losses  \RX\  and  |J?,2|. 
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An  equivalent  expansion  for  ip(r,  z)  is  obtained  by  noting  that 


TABLE  I 

M  ultipath/Angle  Correspondence 


1 

Wo(k j 


l  +  ^\RlR2\<ieiq(M2-2^k‘dz). 


q= i 


(25) 


After  substitution  of  (25)  into  (19),  we  note  that  the  major  con¬ 
tribution  to  the  contour  integral  occurs  within  a  specific  interval 
(fci,fc2)  on  the  real  axis  (see,  for  example,  [13,  ch.  4]).  After 
some  algebraic  manipulation,  the  qth  term  of  the  resulting  se¬ 
ries  expansion  for  tp(r,  z)  is 

I  f  2  H^\kr)u1(z)u2(z)\R1R2\qe^1+^2~2^kzdz)kdk. 

Jki 

(26) 


Note  also  that  the  finite  integration  interval  corresponds  to  the 
use  of  only  N  modes  of  a  possibly  infinite  set  in  the  normal 
mode  sum  and  each  term  corresponds  to  a  specific  propagation 
path  in  the  ocean. 

Separation  of  the  gth  integral  into  four  path  types  is  accom¬ 
plished  by  substituting  the  fundamental  wave  solutions  m  and 
u2  from  (16)  into  (26)  and  noting  that  they  are  linearly  depen¬ 
dent  at  an  eigenvalue.  To  simplify  the  notation,  denote  the  phase 
exponents  as 


Neglecting  amplitude  terms,  then 

ui(z)u2(z0)e~<c 

=  (e~*-  +  R2e^)(e-t°  + 

=  +  R2e-^e^  ](e*«  +  RuT*') 

=  [e-^1  +  R2e-2*ce*rl](e*‘ 

But  R2e~2i^c  =  l/i?i  so  we  have 

Ul(z)u2(z0)e~^c  =  *  +  1  /Il\  —  U\(~'^). 


Path 

Source  Angle 

Receive  Angle 

(1) 

down 

up 

(2) 

up 

up 

(3) 

down 

down 

(4) 

up 

down 

corresponds  to  a  downward  source  angle  and  upward  receiver 
angle,  etc.  The  path  correspondence  is  summarized  in  Table  I . 

Evaluating  (26)  by  the  method  of  stationary  phase  [13],  we 
note  that  for  fixed  receiver  depth  zr,  each  integral  evaluated  at 
range  r  is  of  the  form 

I(r)=  f  2  g(k)e~ifWdk  (28) 

Jki 

where  g(k)  is  slowly  varying  in  comparison  to  e~if(^k\  and  f(k) 
has  continuous  derivatives  to  the  second  order.  Thus  f(k)  can 
be  expanded  about  an  arbitrary  point  k*  within  the  integration 
interval 

m  =  im + f(k*)(k  -  k*) + 1/2  -  k*)2 + . . . 

If  a  stationary  phase  point  k*  with  f(k*)  =  0  can  be  found 
within  the  interval,  (28)  simplifies  to 


I(r )  =  <i(k*)(rif^  [  2  e^2/"0^Hfc-fc*)2d]fc.  (29) 

Jkx 

Assuming  f"(k*)  >  0,  the  change  of  variable  x  = 

y/l/2f"(k*)(k  -  k*)  yields 


I(r) 


g(k*)e^if(k~) 

JW 


(30) 


which  can  beintegrated  numerically.  In  cases  where  the  integra¬ 
tion  limits  can  be  considered  i  nfinite  (i  .e.,  in  intervalscontaining 
a  stationary  phase  point  A*  and  integration  limits  sufficiently  far 
from  k*),  the  integration  can  be  eliminated  by  noting  that 


(31) 


If  f"(k*)  <  0,  the  solution  is  the  complex  conjugate  of  that 
for  f"{k*)  >  0.  In  some  cases  the  solution  is  extrapolated  into 
diffraction  regions,  and  when  there  is  no  stationary  phase  point 
within  the  interval,  integration  by  parts  is  used  to  approximate 
the  basic  integral. 


The  result  of  multiplying  the  two  factors  yields  the  four  prin¬ 
cipal  terms 

e~i^sr  +  R1e-**e-*‘  +  ^  e*' +  e*”.  (27) 

(1)  (2)  (3)1  (4) 

These  four  terms  correspond  to  the  four  possible  combinations 
of  vertical  path  direction  at  the  source  and  receiver,  i.e.,  path  (1) 


IV.  Extension  to  Horizontally  Variable  M  edia 

The  extension  of  the  results  of  Section  III  to  media  in  which 
environmental  properties  may  vary  with  range  is  accomplished 
in  two  stages.  The  first  generalizes  the  calculation  of  the  vertical 
phasefunction  (1)  and  its  derivatives  to  accommodate  horizontal 
changes  in  bottom  slope,  sound  speed,  and  bottom  type  and  the 
second  accounts  for  the  change  in  horizontal  wavenumber  with 
range. 
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A.  Generalization  of  Vertical  Phase  Calculations 

Following  the  notation  of  Section  III  (28),  we  generalize  the 
phase  calculations  to  account  for  horizontal  variations  in  bottom 
depth,  bottom  type,  and  sound  speed.  As  before,  f(k)  corre¬ 
sponds  to  range,  and  f"(k)  to  the  derivative  of  range  with  re¬ 
spect  to  horizontal  wavenumber  A,  which  can  be  approximated 
locally  by  (Ar/Afc). 

In  Section  III  [e.g.,  (25),  (26)],  boundary  losses  were  raised 
to  a  power  q  which  represents  cycle  number.  This  assumes  the 
loss  on  each  boundary  interaction  is  the  same,  and  the  cycles 
of  each  path  are  homogeneous.  When  bottom  and  sound  speed 
parameters  vary  with  range,  this  may  not  be  the  case.  The  no¬ 
tion  of  cycle  carries  over  as  traversal  of  a  wave  from  one  upper 
turning  point  or  reflection  to  the  next,  but  the  paths  traversed 
within  these  cycles  may  vary. 

To  account  for  horizontal  variations  in  sound  speed,  phase 
velocities  are  corrected  according  to  Snell’s  Law.  To  account 
for  variable  bathymetry,  it  suffices  to  detect  when  the  bottom 
is  encountered  within  a  given  rectangle,  and  adjust  the  ray  cal¬ 
culations  to  account  for  the  local  bottom  slope.  On  passing  a 
range  at  which  the  bottom  type  changes,  bottom  loss  calcula¬ 
tions  utilize  the  properties  of  the  new  bottom  type.  As  rays  are 
traced  through  the  medium,  the  values  of  the  phase  derivatives 
(4),  (2),  and  f"  «  Ar/Ak  are  accumulated  as  a  function  of 
range  and  used  in  the  minimum  phase  expansion  of  the  mode 
integrals  (30). 

B.  Adjustment  for  Horizontal  Wavenumber  Differential 

One  final  adjustment  is  required  to  extend  the  mode  integrals 
(30)  into  horizontally  variable  media.  Note  that  when  the  ray 
calculations  are  adjusted  for  SSP  and/or  bottom  depth  changes, 
a  given  phase  velocity,  Cy,  or  equivalently,  horizontal  grazing 
angle,  8,  is  adjusted  in  accordance  with  either  Snell’s  Law 
or  specular  reflection,  respectively.  This  adjustment  results  in 
association  of  the  given  ray  path  with  a  different  horizontal 
wavenumber,  k. 

In  terms  of  the  integral  form  (26),  this  adjustment  necessitates 
a  correction  to  the  integrand  to  account  for  the  change  of  vari¬ 
able  from  the  original  wavenumber  k0  to  the  local  wavenumber 
k(k0).  Considering  the  mode  summation  (24),  which  the  inte¬ 
gral  approximates,  the  change  in  wavenumber  corresponds  to 
a  coupling  of  the  normal  modes  between  the  two  adjacent  ver¬ 
tical  regions.  A  discussion  to  demonstrate  thechangeof  variable 
which  accounts  for  this  coupling  is  as  follows. 

Consider  the  integral  form  (26)  which  corresponds  to  the 
qth  term  in  the  series  expansion  for  <p(r,z),  and  assume  the 
asymptotic  form  of  the  Hankel  function,  and  exponentials  as 
fundamental  solutions  to  the  Green’s  function  (16).  Assuming 
boundary  losses  (15)  are  included  in  the  amplitude  functions 
denoted  by  v(z,k),  and  denoting  the  elapsed  vertical  phase 
along  the  ray  path  r  including  reflection  effects  by  £r(fc),  we 
formulate  an  individual  mode  integral  for  a  receiver  at  range  r 
and  depth  as 

I(r )  =  [  v(z0,  k0)v(zr,  feo)<r<«r(fc°)+fcor)  y/kodk0  (32) 


phase  increment 


The  final  factor  of  in  the  integrand  is  a  combination  of 
k0  and  1/Vko  which  two  factors  appear  at  two  different  points 
in  the  development.  The  first  is  from  the  factor  2k0dk0  which 
resulted  from  a  change  of  variable  from  k%  to  k0  in  considering 
the  A-plane  equivalent  to  the  contour  integral  in  the  k2  plane 
(see  [1,  Section  2.1]),  i.e.,  d(kl)  =  2k0dk0.  The  second  l/v% 
resulted  from  use  of  the  asymptotic  form  of  the  Hankel  func¬ 
tion.  In  passing  to  the  range  dependent  formulation,  the  first  of 
these  factors  remains  the  same,  but  the  second  must  account  for 
variation  in  the  local  wavenumber. 

We  reformulate  the  phase  integral  (3)  by  equating  spatial  and 
temporal  phase 

cut  =  j  kzdz  +  kr  =  f(k)  +  kr  (33) 

and  substitute  £r(fc0)  =  utr(k0)  -  k0rr(k0)  so  that  (28)  be¬ 
comes 

/„— iujtr(ko) 

v(z0,  k0)v(zr ,  ko)  f^-^k0dk0. 

(34) 

To  account  for  the  change  in  horizontal  wavenumber  from 
source  to  receiver,  some  of  the  ko s  must  be  replaced  by  k(k0), 
the  local  wavenumber  at  field  point  (r,zr)  for  the  ray  with 
wavenumber  k0  at  the  source 

/„— iutr(k0) 

A(zo,  zr,  k0)  /—--e-ik^r-rr^k0dko 
vk(ko) 

(35) 

where  A(z0,zr,k0)  =  v(z0,k0)v(zr,k(k0)). 

The  sketch  of  a  sample  path  in  Fig.  1  indicates  the  compo¬ 
nents  of  phase  for  a  receiver  at  range  r.  Clearly,  the  amplitude 
factors  associated  with  source  and  receiver  are  to  be  evaluated 
at  k0  and  k(k0),  respectively,  and  the  sketch  also  indicates  the 
use  of  k0  and  k(k0)  in  the  phase  exponents. 

To  justify  the  value  ^k(k0)  in  the  denominator,  consider 
the  accumulated  transmission  coefficients  across  the  vertical  re¬ 
gions  from  source  to  receiver  at  ranger.  Consider,  for  example, 
a  transition  at  range  r0  from  region  i  to  i  +  1  with  horizontal 
wavenumbers  h  and  ki+1,  respectively,  as  depicted  in  Fig.  2. 
Assuming  an  incident  plane  wave  of  unit  amplitude,  continuity 
of  the  wave  function  and  first  derivative  at  r0  yields 


where  the  subscript  r  denotes  accumulation  along  the  path, 
and  k0  denotes  the  horizontal  component  of  wavenumber  at  the 
source  position. 


eiki( J’-J’o)  _|_  j^e-iki(r-r0)  _  rp  gifc<+i(f— *«) 


r„)  _  ikiRie-ik,(r-r0)  =  ik.+lT.  e*+ i(r-ro) 


and 
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(a) 


region  L+l 


Fig.  2.  Horizontal  transmission  coefficient. 
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1  +  B,j  —  Tj  and 

ki(l  -  Ri)  =  ki+iTi. 


oc 
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CL 


(b) 


(C) 


Fig. 4.  ExampleCaseB:TL  versus  range  for  a  350-H  z  source  i  n  a  convergence 
zone  environment. 


Solving  for  Tit  we  obtain  the  magnitude  of  the  transmission 
coefficent 


T,  = 


2k  i 


ki  +  h+i 


We  perform  some  algebraic  manipulation  to  approximate  T;  in 
a  more  convenient  form 


Ti  = 


2kn 


h+1  +  ki  1  +  'Li±i  2  + 


ki  “  1  ki 

Thus,  to  first  order,  we  can  approximate 

„  1  Hki 


1  + 


2  ki 


Fig.  3.  Example  Case  A:  TL  versus  range  for  a  3520-Hz  source  in  a  surface 
duct  environment. 


1  + 


k 


'i+ 1 


Over  N  vertical  regions,  the  overall  transmission  coefficient  T 
is  approximated  by  forming  the  product 


N 


r=n* 


2=0 


fcjV-1 

k'N 
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Fig.  5. 


Range  (kyd) 


Example  Case  C:  TL  versus  range  for  a  3860-Hz  source  in  shallow  water. 


Noting  that  kN  in  this  case  is  k(k0),  then,  we  see  that  the  am¬ 
plitude  of  the  radial  propagation  is  altered  in  such  away  that  the 
original  v^o  in  the  denominator  is  replaced  by  ^k(k0). 

To  apply  the  method  of  stationary  phase  to  the  range  depen¬ 
dent  form  of  the  mode  integral  (26),  following  the  use  of  the 
Fresnel  integral  form  (28),  we  have 


I(r)  =  A(zq,  zr,  kg)- 


f  dr  /  dk 
dk0  \l  dko 


^-i$(k%,r)  /  e~ix2 


where  .4(20,  zr,  /,:,*)  =  11(20,  kQ)v(zr,  k^),  and  4>(A,o,r)  = 
(Cr(^'o)  +  fc(fc5)r.  Noting  that  ( dk/dko )  =  (tan9r/tan9or): 

where  90r  and  6r  are  the  ray  angles  at  the  source  and  receiver, 
and  using  the  asymptotic  form  for  the  amplitude  functions 


v(zr,kl)  = 


/kz(zr,k%)  tan  0Or 


I{r )  =  v(zQ,k*0) 


ygl  ( 

ltai\6rk(k*0 


-i£(k*,r)  J  e~ix 


TABLE  II 

Bottom  Sediment  Parameters  for  Example  Case  D 


Range 

0 

3 

6 

7.5 

12.09 

Sthiek 

1.091 

1.091 

1.091 

1.091 

1.091 

Xratio 

1.009 

1.0146 

1.0112 

1.0094 

1.0015 

Xgrad 

.634 

.634 

.634 

.634 

.634 

Beta 

25 

25 

25 

25 

25 

Satten 

.06477 

.06477 

.06477 

.06477 

.06477 

Sagrad 

0 

0 

0 

0 

0 

Sdense 

1.8 

1.8 

.8 

1.8 

1.8 

Xdense 

1.8 

1.8 

1.8 

1.8 

1.8 

Xthick 

0 

0 

0 

0 

0 

Baslos 

.069 

.069 

.069 

.069 

.069 

Finally,  recognizing  that  (l/\A’(^o)tan(9,,)  is  asymptotically 
v(zr,  k),  we  have  the  required  result 

I(r )  =  Aiz0iZr,k(ko))V^^^er^^  [  e~ix'2  dx. 

w  —  J 

y  dko 

(36) 

This  result  is  the  range  dependent  analog  of  (30).  The  ampli¬ 
tude  factors  v  are  computed  as  before  but  as  a  function  of  the 
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Fig.  6. 
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Example  Case  D:  TL  versus  range  for  a  142-Hz  source  in  an  environment  with  horizontal ly  variable  bottom  depth  and  sound  speed. 


local  wavenumber  at  source  and  receiver.  The  Fresnel  integral 
is  also  computed  as  in  the  original  derivation,  i.e.,  by  substi¬ 
tution  of  (31)  in  the  case  of  infinite  integration  limits,  or  by  an 
approximation  based  on  integration  by  parts.  The  remaining  fac¬ 
tors  involving  the  (range  dependent)  accumulation  of  the  phase 
functions  al  ong  the  path  r,  and  the  hori  zontal  wavenumber  asso¬ 
ciated  with  the  stationary  phase  point  comprise  the  changes 
required  to  extend  the  solution  to  horizontally  variable  media. 

V.  Discussion  of  Results 

Thefol  lowing  provides  samples  of  calculations  made  with  the 
integrated  M  ode  model  described  herein,  with  comparisons  to 
predictions  computed  using  other  models  and  measured  results. 
In  each  case,  measured  results  or  other  model  calculations  are 
shown  as  dots,  and  the  prediction  of  the  I  ntegrated  M  ode  model 
is  represented  by  a  solid  line. 

A.  Surface  Duct  Propagation 

The  first  sample  case  involves  propagation  of  a  3520-Hz 
signal  in  a  surface  duct  [14],  The  field  is  dominated  by  a  small 
number  of  low  order  modes  and  both  measurements  and  model 
calculations  are  restricted  to  propagation  within  the  duct.  The 


TABLE  III 

Bottom  Sediment  Parameters  for  Example  Case  E 


Sthick 

1.32  sec 

Xratio 

0.999 

Xgrad 

1.3 

Beta 

-0.5 

Satten 

0 

Sagrad 

0.0001 

Sdense 

1.44 

Xdense 

1.63 

Xthick 

0.65 

Baslos 

0.45 

source  and  receiver  are  in  the  duct  at  depths  of  20  and  50  feet, 
respectively,  and  the  wind  speed  is  5  knots.  The  upper  1000 
feet  of  the  SSP  and  the  results  for  Sample  Case  A  are  shown  in 
Fig.  3(a)  and  (b),  respectively. 

B.  Deep  Water  Convergence  Zone  Propagation 

The  second  sample  case  involves  propagation  of  a  350-Fiz 
signal  in  a  convergence  zone  environment  [15].  The  source,  re¬ 
ceiver,  and  bottom  depths  are  400,  600,  and  16  920  feet,  respec¬ 
tively,  and  the  wind  speed  is  4  knots.  The  bottom  was  mod- 
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(a) 


Sound  Speed  (ft/s) 
(b) 


Fig.  7.  Sample  Case  E:  TL  versus  range  for  a  3500-Hz  source  in  the  Norwegian  Sea. 


eled  using  a  plane  wave  reflection  coefficient  with  the  loss  per 
bounce  as  shown  with  the  SSP  and  the  results,  in  Fig.  4(a),  (b), 
and  (c),  respectively. 

C.  Shallow  Water  Propagation 

The  third  sample  case  represents  propagation  of  a  3860-Hz 
signal  in  a  shallow  water  environment  with  a  single  downward- 
refracting  SSP  and  highly  variable  downslope  bathymetry  [16], 
The  source  and  receiver  were  both  at  a  depth  of  25  feet.  The 
bottom  was  modeled  using  a  plane  wave  reflection  coefficient 
with  the  loss  per  bounce  as  shown  in  Fig.  5(a).  The  bathymetry 
profile  for  Sample  case  C  and  an  overlay  of  the  sound  speed 
variation  with  depth  are  shown  in  Fig.  5(b),  and  the  results  for 
Sample  Case  C  are  presented  in  Fig.  5(c). 

D.  Propagation  in  an  Environment  With  Plorizontally  Variable 
Sound  Speed,  Bottom  Depth,  and  Bottom  Type 

The  fourth  sample  case  represents  propagation  of  a  142-Hz 
signal  in  a  shallow  water  environment  with  horizontally  vari¬ 
able  sound  speed,  bottom  type,  and  bathymetry  [17],  The  SSPs 
at  ranges  of  0,  3,  6,  7.5,  and  12.09  km  are  plotted  with  the  ups- 
lope  bathymetry  forSamplecaseD  in  Fig.  6(a).  The  bottom  sed¬ 
iment  parameters  by  range  are  given  in  Table  II,  where  Sthick  is 


sediment  thickness  represented  by  two-way  travel  time  through 
the  sediment  in  seconds,  X  ratio  is  the  ratio  of  sediment  to  water 
column  sound  speed  at  the  sediment  interface,  Xgrad  is  the 
initial  sediment  sound  speed  gradient  in  1/s,  Beta  is  a  sedi¬ 
ment  SSP  curvature  parameter,  Satten  is  initial  sediment  atten¬ 
uation  in  dB/m/kFIz,  Sagrad  is  sediment  attenuation  gradient  in 
dB/m/kFIz/m,  Sdense  is  sediment  density  in  g/cc,  Xdense  and 
X  thick  are  the  density  and  thickness  of  a  hypothetical  thin  layer 
at  the  sediment/water  column  interface  in  g/cc  and  m,  respec¬ 
tively,  and  Baslos  is  the  basement  reflection  coefficient.  The 
source  and  receiver  are  both  at  a  depth  of  18  meters. 

The  results  for  Sample  Case  D  are  presented  in  Fig.  6(b). 

E.  Comparison  to  Navy  Standard  Grab  Model 

The  fifth  sample  case  represents  propagation  of  a  3500-FH z 
signal  in  the  Norwegian  Sea  in  winter  with  results  compared 
to  those  of  the  Navy  Standard  Gaussian  RAy  Bundle  (GRAB) 
model  of  Weinberg  [5].  The  bottom  depth  is  10  500  feet  and  the 
source  and  receiver  are  at  depths  of  325  and  850  feet,  respec¬ 
tively.  Bottom  parameters  are  given  in  Table  III  with  parameters 
and  units  as  defined  in  Sample  Case  D,  and  the  wind  speed  is 
16  knots.  The  SSP  and  results  are  shown  in  Fig.  7(a)  and  (b), 
respectively. 
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VI.  Summary 

The  multipath  expansion  method  of  solving  the  Helmholtz 
wave  equation  to  describe  the  underwater  sound  field  for  a  fixed 
point  source  in  a  plane  multilayered  medium  has  been  discussed. 
This  approach  has  been  extended  to  accountfor  horizontal  vari¬ 
ations  in  bottom  depth,  bottom  type,  and  sound  speed  in  the  sta¬ 
tionary  phase  approximation.  A  comparison  of  calculations  with 
both  measured  data  and  results  of  the  N  avy  standard  model  cov- 
ering  a  range  of  frequencies  and  environmental  variations  has 
been  presented. 
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